Thermodynamics of surface tension: application to electrolyte solutions 
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In this contribution to the special issue of the Journal of Statistical Physics dedicated to Michael 
Fisher on his 70'th birthday, I shall review two thermodynamically distinct routes for obtaining the 
interfacial tension of liquid-vapor interfaces in mixtures. A specific application to the calculation of 
excess surface tension of aqueous electrolyte solutions will be presented. 
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I. INTRODUCTION 

It is a pleasure to dedicate this review to Michael 
Fisher on the occasion of his 70th birthday. The paper 
deals with two subjects which are dear to Michael's heart, 
the Thermodynamics and the Coulomb Systems. Here I 
will review two thermodynamically distinct routes to sur- 
face tension of aqueous electrolyte solutions. The first, 
grand canonical route, goes all the way to the pioneer- 
ing work of Gibbs on the foundations of thermodynamics 
and statistical mechanicsEI and is the usual method used 
by physical chemists. Application of the Gibbs adsorp- 
tion isotherm to the calculation of surface tension otelec- 
trolyte solutions started with the works of Wagnero and 
Onsager and Samarasa early in the 20th century. The 
second, canonical method, has been introduced recently 
and waa.£pund to work very well for symmetric 1:1 elec- 
trolytes! 



II. GRAND-CANONICAL ROUTE 

Many elementary thermodynamics and statistical me- 
chanics texts neglect to deal with all the subtleties lead- 
ing to the Gibbs adsorption isotherm. While these do not 
seriously affect the results for surface tension of surfac- 
tant solutions, the increase of interfacial tension of water 
due to salts is proportionately much smaller. Thus, a 
particular care must be taken with the thermodynamics 
in order to account for all the relevant contributions. We 
start, therefore, by reviewing the thermodynamics lead- 
ing to the Gibbs adsorption isotherrnErQ. 

Consider an aqueous solution in equilibrium with va- 
por, Fig 1. The bulk density of liquid is p[ and the bulk 
density of vapor is p\. The variation of density is con- 
fined to the interfacial region a located between z a and 
Zb, where z is the axis perpendicular to the interface. 
The thermodynamic equilibrium requires the constancy 
of chemical potentials of solute and solvent and of the 
pressure throughout the system. The differential inter- 
nal energy of the interfacial region a is 




interfacial region o" 



FIG. 1: Schematic density profiles inside a solution: pi(z) is 
the characteristic density variation of water across the liquid- 
vapor interface, pi(z) is the concentration profile of solute. 
The hypothetical dividing surface separating the liquid from 
the vapor is located at z = Za- The discontinuous "bulk" pro- 
files for solvent and solute are pi{z) and p2(z), respectively. 
Note that differently from non- ionic solutes for which pi{z) 
extends all the way into the vapor phase, because of a fa- 
vorable gain in solvation free energy aqueous electrolytes are 
completely confined to the liquid phase. 



where T, P, S,V, A, 7, pn, Ni are the temperature, pres- 
sure, entropy, volume, area, surface tension, chemical 
potential, and the number of particles of type i. The su- 
perscript a stands for the interfacial properties. The sum 
runs over all the species, solute and solvent. Since the in- 
ternal energy E is an extensive function of {S,V, A, Ni}, 
application of Eulcr theorem for first order homogeneous 
functions allows the integration of Eq. (m yielding, 



TS a - PV a + jA a + mN° 



(2) 



As usual the Gibbs free energy for the interface is defined 
through the Legendre transform 



dE° 



TdS a - PdV a + -id A" + 22 frdNf , (1) 



G° 



TS a + PV° ~ j A 



(3) 
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which after substitution of internal energy, Eq. 
duces to 



(4) 



The Gibbs free energy is a natural function of 
{T, P, 7, Nf } and its differential is 

dG° = -S a dT + V a dP- A a d 1 + ^^dNf . (5) 

On the other hand differentiating Eq. (Q), we find 

dG a = MdN° + J2 N i d »i ■ ( 6 ) 
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Comparing Eq. (||) with Eq. (^) we are led to a Gibbs- 
Duhcm-likc equation for the interface, 



S a dT - V a dP + A a d 1 + ^ Nf dp, 



(7) 



Now lets consider in more detail a two component sys- 
tem. In this case Eq. (0) simplifies to 



S°dT - V°dP + A a d-f + N[dm + N%dfi 2 







(8) 



The chemical potentials are uniform throughout the sys- 
tem and their variations can be obtained by considering 
the bulk liquid phase. Remembering that a change in 
concentration of solute affects the chemical potentials of 
both solute and solvent, we have 



d/j,i = —s\dT + v\dP ■ 



and 



dfi 2 



-S2dT + v 2 dP 
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T,P 



dc b 



(9) 



(10) 



T,P 



where the si and v\ are the partial entropy and volume 
per particle of solvent, s 2 and v 2 are the partial entropy 
and volume per particle of solute, and c& is the concentra- 
tion of solute, all the values taken inside the bulk liquid 
phase. Substituting Eqs. (Oft and (|l(]) into Eq. (ph, the 
Gibbs-Duhem equation for the interface becomes, 



(S a - N[si - N^s 2 )dT - 
{V a - Nfvi - N£v 2 )dP + 

'd^A + 

dc b J T p 



(11) 



A°dj 



N? 



dfi 2 
dc b 



T,P 



dch = 



For a truly two component system pressure cannot be 
held constant if T or Cb is varied. However, presence of 
an inert gas does not significantly affect aqueous,surface 
tension, so that in practice P can be kept fixedcl. With 
the help of the Gibbs-Duhem equation for the bulk liquid, 



Ni 



dc b 



N 2 



T,P 



djJ-2 
dc b 



= 0, 



(12) 



where Ni and N 2 are the bulk numbers of solvent and 
solute molecules, Eq. (If) reduces to two equations, 



df 



Cb,P 



(13) 



and, 



d-f 

9(^2 / T,P 



Ni 



(14) 



T,P 



where Ti = Nf j A a and T 2 = Nf/A a . Equation Ql£ 
states that the variation of surface tension with respect 
to temperature is minus the excess entropy, compared 
to the entropy content of the same amount of material 
inside the liquid phase. The equation (Q), called the 
Gibbs adsorption isotherm, shows that the change of sur- 
face tension with respect to chemical potential of solute 
is minus the excess of solute inside a, over the amount 
which would be present for the same quantity of solvent 
in a bulk liquid phase. Note that the term "excess" is 
used both when the interfacial region has higher (posi- 
tive excess) or lower (negative excess) concentration of 
solute, as compared to the bulk liquid phase. Since the 
thermodynamic stability requires dfi 2 /dcb > at fixed 
temperature, Eq. (Q) shows that a positive surface ex- 
cess of solute leads to a lower interfacial tension, while a 
negative surface excess of solute results in a higher sur- 
face tension. Finally, we note that the Eq. (|l4j) is invari- 
ant with respect to the specific location of the dividing 
surfaces z a and long as they completely enclose the 

region of strong density variation. 

The surface tension can be obtained by integrating 
Eqs. ( |T3| ) and (|T4|). Of course, this requires a specific 
microscopic model which would allow us to calculate the 
excess entropy and the excess amount of solute inside a. 
Such modaLcan, in principle, be provided by statistical 
mechanicsau. 

In the language of statistical mechanics the calcula- 
tions based on Eqs.([l3]) and ( |T4| ) are intrinsically grand- 
canonical. The number of particles inside a fluctuates 
and is determined by the condition of equilibrium be- 
tween the interfacial region and the bulk phases. Unfor- 
tunately to actually solve a statistical mechanical model, 
one is invariably forced to make approximations. It might 
be, therefore, worthwhile to explore different routes to 
surface tension, each one requiring different approxima- 
tions. Of course, in an exact calculation all the ther- 
modynamic routes will lead to the same result. With 
approximate theories this, however, is no longer the case 
and some routes can be significantly better than others. 
With this in sight_|We have explored the canonical route 
to surface tensionQB. 
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III. CANONICAL ROUTE 

Suppose that the system depicted in Fig. 1 which con- 
tains 



Ni = A 



particles of solvent and 



N, = A 



H 



pi(z)dz 



H 



p 2 (z)dz 



(15) 



(16) 



particles of solute is confined to a cylindrical box of height 
H and cross-sectional area A. We shall define the "bulk" 
density of solvent and solute inside the liquid phase as 
Pi = Pi{0) and p\ = ^2(0), respectively; and the "bulk" 
density of solvent and solute inside the vapor phase as 
Pi = Pi (H) and p 2 — p 2 (H) . The internal energy of the 
whole system is 



E = TS - PV + 7 A + fiiNi + p 2 N 2 



(17) 



Now we would like to ask which part of this energy is 
due to the interface? That is, if instead of a continuous 
density profiles pi{z) and P2{z), we would have two dis- 
continuous "bulk" profiles pf(z) and p 2 (z), depicted in 
Fig. 1, what would be the change in internal energy of 
the system? Although easy to pose, this question carries 
some subtlety. Specifically where should we put the di- 
viding surface? Also, what is the meaning of a discontin- 
uous profile? This latter question is fairly easy to answer. 
To obtain a discontinuous profile we divide the total vol- 
ume into two sub-cylinders of heights z d and H — z d . In 
the first sub-cylinder of height z d and volume V 1 = Azd, 
there will be N\(z d ) = p\V molecules of solvent and 
N 2 (zd) = p l 2 V l molecules of solute. In the second sub- 
cylinder of height H — Zd and volume V v = A(H — z d ), 
there will be N^(zd) — p\V v molecules of solvent and 
N%{zd) = p 2 V v molecules of solute. To keep the uniform 
density distribution characteristic of the bulk phases, we 
impose periodic boundary conditions in the z direction 
on each sub-cylinder. The internal energy of the first 
sub-cylinder, corresponding to the bulk liquid, is then 



Helmholtz free energy is given in terms of the Legendre 
transform of the internal energy, F s — E S ~TS S , Eq. (jgtj ) 
can be rewritten as, 



F s {z d ) - pxNl{z d ) - p2N§(z d ) = ~/A 



(21) 



The statistical mechanics allows us, in principle, to cal- 
culate the excess surface Helmholtz free energy as well as 
the chemical potentials and the surface excess of solvent 
and solute. The thermodynamics relates these to sur- 
face tension through Eq. (pl|). Note that while the three 
terms on the left hand side of Eq. ( |2l| ) are dependent on 
the location of the dividing surface, the right hand side 
does not. Thus, we can fix the position of Zd so that the 
surface excess of solvent is zero, Nl(z d ) = 0. This spe- 
cific choice corresponds to the, so called, Gibbs dividing 
surface. With the location of z^ specified, the surface 
tension becomes 



F°(zC)-p 2 Nj(zC) 
A 



(22) 



Now, lets define the liquid surface excess of solute as 



A l =A [ P2 (z) - p\\dz 
Jo 

and the vapor surface excess of solute as 

r-H 



A V =A[ [p 2 (z) - pl\dz . 



Noting that N 2 = A 1 + A v and 



M2 



dF Lik 
dNi 



T,V> 



dF bulk 



T,V V 



(23) 



(24) 



(25) 



Eq. (p2| ) can be rewritten as 



7 = _ [F(N U N 2 ) - FL lk (N[,N l 2 ) - FZ ulk (N?,N%)- 



dNi, 



8N$ 



(26) 



E l {z d ) = TS l {z d ) - PV 1 + piN[{z d ) + p 2 N l 2 (z d ) , (18) 

and the internal energy of the second sub-cylinder, cor- 
responding to the bulk vapor, is 

E v (z d )=TS v (z d )-PV v +n 1 N?(z d ) + n 2 N$(z d ) . (19) 

Because of the arbitrariness in the location of the dividing 
surface, in general, Ni ^ N{(z d ) + N^(zd) and N 2 ^ 
N 2 (zd) + N 2 (z d ). Therefore the surface internal energy, 
E s (z d ) =E-E l (z d )-E v (zd), is 

E s (z d ) = TS s (z d ) +-/A + piNt(zd) + p 2 m(z d ) , (20) 

where S s {z d ) = S - S l {z d ) - S v (z d ) , Nf = Ax - N[ (z d ) - 
Nf(zd), and iVJ = N 2 - N l 2 {z d ) - N%{z d ). Since the 



where F(Ni, N 2 ) is the total Helmholtz free of the system 
with interface. 

In the thermodynamic limit A 1 /N l 2 < 1 and A v /N% < 
1, and Eq.(]26|) simplifies to 

7 = j[F(N 1: N 2 ) - F l bulk {NlN l 2 + A 1 )- (27) 
i^ fc (AM + A")]. 

This is the principal thermodynamic result of this pa- 
per, which will serve as the starting point for our analy- 
sis of the interfacial tension of electrolyte solutions. The 



fact that Ni = N[ + N% and N 2 



N l 2 



A 1 + A v 

signifies that the calculations based on Eq. ( |27| ) must 
be performed with a fixed number of solvent and solute 
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molecules. This accounts for the adjective "canonical" in 
the name of this thermodynamic route to surface tension. 

For strong electrolytes, we have an additional simplifi- 
cation. Because of the high dielectric constant of water, 
ions inside the liquid have a significantly lower electro- 
static free energy compared to ions inside the vapor. The 
probability of finding an ion inside bulk liquid pi, com- 
pared to its probability of being inside bulk vapor p v is 



El 

Pv 



-0W 



(28) 



where (3 = l/k B T and W is the solvation energy, which 
can be estimated from the Born equation, 



W = 



i 



(29) 



In this formula q is the ionic charge, d is the ionic di- 
ameter, e; is the dielectric constant of liquid, and e v is 
the dielectric constant of vapor. For water ei/e v w 80, 
so that PW can be well approximated by —\ B ei/de v , 
where X B = (3q 2 /ei is the Bjerrum length. For water 
at room temperature and monovalent ions X B = 7.2 A, 
while the characteristic diameter of a hydratcd ion is ap- 
proximately d fts 4 A. Therefore, the concentration of ions 
in vapor can be estimated to be p\ = p\ exp(— 140), which 
is immeasurably low. Furthermore, a rapid decrease of 
the dielectric constant for z > z^ should strongly inhibit 
presence of ions in the interfacial region, confining them 
to z < zf, so that N$ = and A" = 0. 



IV. A SIMPLE MODEL OF AN ELECTROLYTE 

Let us now consider a simple model of an aqueous elec- 
trolyteO confined to a cylinder of cross-sectional area A 
and height H . The N ions will be idealized as hard 
spheres of diameters d, carrying charge +q or — q at their 
center. Suppose that the Gibbs dividing surface is lo- 
cated at the top of the cylinder. To simplify the model, 
we shall further assume that on crossing the Gibbs di- 
viding surface from bellow the dielectric constant drops 
discontinuously from e = 80, characteristic of bulk wa- 
ter, to e = 1, characteristic of vacuum. From the previous 
discussion, the increase in surface tension of water due to 
electrolyte is 
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lim 

A,H,N- 



1 

A 



( pel pel \ 



(30) 



where the thermodynamic limit is taken in such a way 
as to preserve the bulk concentration, Cb = N/2AH. In 
Eq. (|3C|), F ex is the excess Helmholtz free energy due 
to electrolyte in the presence of liquid-vapor interface, 
while F b ex lk is the excess free energy of the same amount 
of electrolyte but dissolved within the bulk liquid phase of 
the same volume. The F b ex [k can be calculated using the 
periodic boundary conditions at the top and the bottom 
of the cylinder. 



As was already mentioned, the decrease in free energy 
due to solvation inhibits entrance of ions into the inter- 
facial region, producing an ion-free layer adjacent to the 
Gibbs dividing surface. The different hydration charac- 
teristics of anions and cations can make their exclusion 
layers to have different width. Let us call the width of a 
cation exclusion layer <5 + and the width of an anion exclu- 
sion layer (5_. Furthermore, lets for the moment neglect 
all the electrostatic and hardcore interactions, beyond 
their contribution to the formation of ion-free layers. In 
this case the total free energy of an electrolyte solution 
inside the cylinder takes a particularly simple entropic 
form, 

F ex = k B TN+ [ln(c+A 3 ) - l] + k B TN_ [ln(c_A 3 ) - l] , 

(31) 

where N+ = N_ = N/2, c+ = N + /A(H - 6+), 
c_ = N-/A(H — 8-), and A is the thermal de Broglie 
wavelength. On the other hand the bulk free energy, 
-^buifc! i s gi ven by exactly the same expression, but with 
c + = c_ = Cb = N/2AH, since the periodic boundary 
condition destroys the ion free layers. Substituting into 
Eq. (ppj) we find a simple result, 



Y x = k B Tc b (5+ + S-) 



(32) 



Neglect of direct electrostatic interactions results in ex- 
cess surface tension, scaling linearly with the concentra- 
tion of elect rolyteBQ. Eq. (32) is quite different from the 



Onsager- Samaras (OS) limiting la- 



ILL 



70 [- In(/sA B /2) - 2 lE + 3/2] 



(33) 



where 7_g — 0.577215665... is the Euler's constant, 70 = 
q 2 Cb/2e, and k — ^/8irq 2 Cb/ek B T is the inverse Debye 
length. The OS limiting law is believed to be univer- 
sally valid for all electrolytes at infinite dilution. In this 
respect, it is supposed to be analogous to the Debye lim- 
iting laws for bulk electrolytes. Comparing Eq. ( |33| ) with 
Eq. ( |32] ) we see, however, that only the first term inside 
the square brackets is universal, while the existence of ion 
free layers leads to corrections which scale linearly with 
the electrolyte concentration. Although, Onsager and 
Samaras tried to extend their theory to finite densities, 
they have overlooked existence of ion free layers and thus 
missed an important contribution to the surface tension. 

The experimental measurements show a purely linear 
dependence of excess surface tension on concentrationEl, 
with the validity of OS limiting law restricted to very 
low electrolyte densities. Neglect of electrostatic inter- 
action is clearly a strong oversimplification. The advan- 
tage of the canonical formalism, besides its thermody- 
namic simplicity, is that the electrostatics can be eas- 
ily rtf^f*! m f° account. Thus, the Debye-Hiickel the- 
orylallJ'LlI for bulk electrolytes can be jegtended to take 
into account the liquid-vapor interfacaju. This calcula- 
tion, which is in excellent agreement with experiments, 
shows that for NaCl, existence of an ion free layer of 
width S + = 6- = 2.125 A, characteristic of ionic hydra- 
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tion radius, accounts for about 30% of the total contribu- 
tion to 7 ea: , for 0.5M to IM solutions. Electrostatics be- 
comes even more important at lower concentrations. For 
solutions with q, < 0.3M the electrostatics is over 80% 
dominant, but the OS limiting law fails to be a reason- 
able approximation until q, < 0.1 M. Furthermore, even 
an extended grand-canonical calculation^ which goes be- 
yond the limiting laws, significantly underestimates the 
role of electrostaticstm. This suggests that the canonical 
route to surface tension of aqueous electrolytes is more 
reliable. 

It would be interesting to see if the electrostatic cal- 
culations for symmetric electrolytes can be extended to 



the asymmetric systems with distinct <5 + and 6- , perhaps 
along the same lines as used recently by Michael Fisher 
and collaborators for the asymmetric bulk electrolyteslia. 



V. ACKNOWLEDGEMENT 

It is a pleasure to acknowledge interesting conversa- 
tions with A. Parsegian. This work was supported in 
part by Conselho Nacional de Desenvolvimento Cienti'fico 
e Tecnologico (CNPq). 



1 J. W. Gibbs, Collected Works, Longmans, Green and Co., 
NY (1928). 

2 C. Wagner, Phys. Z. 25, 474 (1924). 

3 L. Onsager and N. N. T. Samaras, J. Chem. Phys. 2, 528 
(1934). 

4 Y. Levin, J. Chem. Phys. 113, 9722 (2000). 

5 Y. Levin and J.E. Flores-Mena, Europhys. Lett. 56, 187 
(2001). 

6 J.E.B. Randies, in "Advances in Electrochemistry and 
Electrochemical Engineering" (P. Delahay and C.W. To- 
bis, Eds.), Vol.3, p.l, Interscience, NY, 1963. 



7 E. Schmutzer, Z. phys. Chem. (Leipzig) 204, 131 (1955). 

8 N. Matubayasi, H. Matsuo, K. Yamamoto, S. Yamaguchi, 
and A. Matuzawa, J. Colloid Interface Sci. 209, 398 (1998). 

9 P. W. Debye and E. Hiickel, Phys. Z. 24, 185 (1923). 

M.E. Fisher and Y. Levin, Phys. Rev. Lett. 71, 3826 
(1993). 

1 Y. Levin and M.E. Fisher, Physica A 225, 164 (1996). 

2 D.M. Zuckerman, M.E. Fisher, and S. Bekiranov, Phys. 
Rev. E 64, 011206 (2001). 



